/* Phillips curve estimation program for Rotemberg pricing */
/* By Toshihiro OKADA */


//***** Choosing a set of parameters (posterior modes): Basic NK, Standard NK, and EPD-NK models *****//  
@#define data_techzero=1 
//* 0: Basic or Standard NK                      *//
//* 1: EPD-NK                                    *// 

//********** Model choice ********* //
@#define tech =1
//* 1: EPD NK                *// 
//* 0: Basic or Standard NK  *//

//**** Shutting down endogenous technology for EPD-NK model ***// 
@#define ashock=1
//* 1: with endogenous tech.     *//
//* 0: without endogenous tech.  *//
//* 2: arbitrary value for "alfa" (R&D elasticity)  *//

//*** Rotemberg VS Calvo  ***//
@#define roten =1
//* 0: Calve pricing          *//
//* 1: Rotemberg pricing      *//

//***** Price index *****//
@#define pindex=0
//* 0: no price index  *//
//* 1: with price index *//


//***** Params setting for couterfactural analyses ******// 
alfa_n = 150; /* to shut down technology effect */
mode_adj =1.0; /* to choose an arbitrary value of "alfa": only go with ashock==2 */
rho_rot =0.7218; /* to set the value of the Rotemberg pricing parameter value, "rot", to match with the Calvo parameter, "rho", (0.7218: posterior mode of "rho") */ 


//*** Time-to-innovate specification *****//
@#define TI=1
//* 1: Time to Innovate (8 periods)   *//
//* 0: No Time to Innovate (1 period) *//
//********************************//


//******** Conditional varinace; on-off ****//
@#define cndvar=0
/** cndvar=1: for (theoretical) conditional variance calculation **/  
/** cndvar=0: for otherwisee **/
//*******************************************//

//******* Mone Carlo or Theoretical IRFs  ******// 
@#define fig =0
/** fig=1: for theoretical simulation with IRFs&conditional variance **/
/** fig=0: for monte carlo simulations **/  
//******************************************************************************//

//@# define std_spec=1


//// Line specification for figures ////////////
lsp =2; /* setting a line specification for figures: 1='-', 2='--', 3=':' and 4='-.' */



//***********************************************//
//******** Loading data ************************//
//**********************************************//
@#if data_techzero==0

//* Without Price Indexation (Basic NK) *// 
@#if pindex==0
load('subdraws_BNKphi60stv90m3_mode.mat', 'aa0');
@#endif

//* With Price Indexation (Standard NK) *// 
@#if pindex==1
load('subdraws_BNKphi60stv90index_mode.mat', 'aa0');
@#endif

@#endif

@#if data_techzero==1

load('subdraws_phi60stv90mm5_mode.mat', 'aa0');

@#endif
//***********************************************//

//****************************************************//
//**** Estimation sample size: large and small ******//
//**************************************************//
@#define large_sample=0
//** For the large sample case, choose large_sample=1 and for the small sample case, choose large_sample=0 **//
sim_period=190; //** Set this to 7000 for the large sample case and 190 for the small sample case **//
cut_period=400; //** DO NOT CHANGE **//



//** (Not Used )  **//
 @#define merror=7


//******* IV specification & its Standard errors calculation method ***********//
@#define just_id=0
//* just_id=1: just identified IV, just_id=0: over identified IV calculation +//
@#define nw_stdr=1
//* nw_stdr=1: Newey_West std errors , nw=stdr=0: normal std errors *//
L=12; 
//* L = # of lags for of Neway=West standard errors: See Gali&Gartler (JME1999) *//

//*** With/Without constant terms for IV *****************//
@#define nocon=1
//** nocon=0:with constant, nocon1=1:no constant *//

//**** NOT USED: Scale adjubstment for errors **///
@#define scl_adj=0


//********** Measure ment erros in marginal costs: 1: NO NEED TO CHANGE *************************************//
//**** Correlation adjubstment for errors **///
 @#define corr_spec_mc=1
 @#define corr_spec_fpai=1 
//** 1 if randwithcorr.m (corr & standard err rest) is used and "1" should be chosen even with zero mesurement error case  **//


//********** Measure ment erros in marginal costs: 2: NEED TO CHANGE IN THE CASE OF MEASUREMENT ERRORS *************************************//
crr_mc= 1.000; /*�@set it to 0.9 for measument errors in MC. For measument errors in mc (correaltion with true MC): corr(mc,mce):correlation coefficient between mc and (mc+mee) */ 
crr_fpai = 1; /*�@Not used */
mc_adj=1.0; /*�@Not used */
mc_adj2=1.0; /* set it to 0.9 for measument errors in MC. corr(fmce, fmc) where fmc is expected mc and fmce is fmc with errors*/


//**** MC adjustment as an instrument: NO NEED TO CHANGE *****///
@#define mc_adj=2
//** 0 = no adjustment, 1=mce, 2=fmc ****************//


//*** Extra lags of ppai (# of lags of inlation) ***
@#define ext_lg=0

adj_=0;/*�@Not used */ 
pxcrr_=1;/*�@Not used */ 



@#if tech==0
stde_mc = 0.0309; /*�@Not used */
stde_fmc= 0.0228; /*�@Not used */
stde_fpai = 0.0077; /*�@Not used */
stde_pai = 0.0094; /*�@Not used */

@#endif
@#if tech==1
stde_mc =   0.0462; /*�@Not used */
stde_fmc= 0.0337; /*�@Not used */
stde_fpai =  0.0139;  /*�@Not used */
stde_pai =  0.0159; /*�@Not used */

@#endif


//********************** NOT USED *******************************//
mee_std = 0; 
meee_std =0;   
me_std=0; 
scl1= ( (stde_mc)^2/(stde_mc^2 + mee_std^2) )^(1/2); 
scl2= ( (stde_fpai)^2/(stde_fpai^2 + meee_std^2) )^(1/2); 
//**********************************************************//



///**************** Number of simulation replications  **********************//
iia=1; /* Start of random shock seed # */
iib=10000; /* 10000. End of random shock seed # */
ja=1; /* Not USED */
jb=1; /* Not USED */



//** Shock Specification for Monte Carlo **//
@#define ushock=1 
@#define news=1  
@#define dshock=0 
@#define cshock=1 
@#define hshock=1 
@#define ishock=1 
@#define fshock=1 
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//* ushock: ushock=0 shut down technology(u)shocks, 2: news: news=0, shut down  news shock but leave u shocks, cshock: consumption shock, hshock: hour shocks, ishock:interest rate shock *// 
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


@#define eshock=16  
/* NOT USED*/


@#define fullshock=1 
/* NOT USED*/
////////////////////////////
//* 1: full shocks   // 
//////////////////////////


@#if fullshock==1
 @#define shki =1
 @#define shku =1
 @#define shkuc =1
 @#define shkuh =1
 @#define shkue =1
 @#define shkue1 =1
 @#define shkue2 =1
 @#define shkue3 =1
 @#define shkue4 =1
 @#define shkue5 =1
 @#define shkue6 =1
 @#define shkue7 =1
 @#define shkue8 =1
 @#define shkue9 =1 
  @#if eshock==16
   @#define shkue10 =1 
   @#define shkue11 =1 
   @#define shkue12 =1 
   @#define shkue13 =1 
   @#define shkue14 =1 
   @#define shkue15 =1 
  @#endif 
  @#if eshock==20
   @#define shkue10 =1 
   @#define shkue11 =1 
   @#define shkue12 =1 
   @#define shkue13 =1 
   @#define shkue14 =1 
   @#define shkue15 =1
   @#define shkue16 =1 
   @#define shkue17 =1 
   @#define shkue18 =1 
   @#define shkue19 =1 
  @#endif
  @#define shkuc =1
  @#define shkuh =1 
  @#define shkf =1 
 @#endif 



for j=ja:jb

@#if tech==1
var fp, ffp, ffpai,pffpai, y_g, c_g, tiv_g, h_g, w_g,  gg, ggn, uf, pai,y,c,k,iv,h,w,r,p,ppai, ppai2,ppai3, ppai4,m,pf,i,lamda,u,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,pfn,in,lamdan,yg,hg,fpai,mcgap,a,rd,rdn,an,nf,nfn,pinf,uc,um,uh,tiv, tivn ynobs,mcgape, mc, mce,fmc,pfmc, fpaie, pfpaie, paie, paiee,ppaie, inftr, tfp, tfpu, tfpn, errmc, shk_me,shk_mee,shk_meee;
@#if eshock==16
varexo ei eu euc euh ef eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue14 eue15 eue me mee meee;
@#endif
 @#if roten==0 
 parameters  tau, pxcrr, adj, hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, psi, epsilon,   alfa, rhoum,rhouc, rhouh, rhou,rhoud, rhouf, sigc, sigh, sigm, ss, rho_eta, eta1, eta2, eta3, eta4,eta5,eta6,eta7,eta8,iss, rss, yss, kss, css,wss,mss, rdss, pfss,  hss, ivss, ggss, lamdass, nfss,x;
 @#endif
 @#if roten==1 
 parameters  tau, pxcrr, adj, hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rot, psi, epsilon,   alfa, rhoum,rhouc, rhouh, rhou,rhoud, rhouf, sigc, sigh, sigm, ss, rho_eta, eta1, eta2, eta3, eta4,eta5,eta6,eta7,eta8,iss, rss, yss, kss, css,wss,mss, rdss, pfss,  hss, ivss, ggss, lamdass, nfss,x;
 @#endif

@#endif

@#if tech==0
var fp, ffp, ffpai, pffpai,y_g, c_g, iv_g, h_g, w_g,  gg, ggn, uf, pai,y,c,k,iv,h,w,r,p,ppai, ppai2,ppai3, ppai4,m,i,lamda,u,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,in,lamdan,yg,hg,mcgap,uc,um,uh, ynobs, fpai,tq,mcgape, mc, mce,fmc, pfmc,fpaie, pfpaie, paie, paiee,ppaie, tfp, tfpu, tfpn errmc, shk_me,shk_mee,shk_meee;

@#if eshock==16
varexo ei eu euc euh ef eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue14 eue15 eue me mee meee;
@#endif
//parameters tau, pxcrr, adj, hb , xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, rhoum,rhouc,rhouh, rhou,rhoud, rhouf, sigc, sigh, sigm, ss, iss_noa, rss_noa, yss_noa, kss_noa, css_noa,wss_noa,mss_noa,  hss_noa, ivss_noa, ggss_noa, lamdass_noa ;

 @#if pindex==0
  @#if roten==0 
  parameters tau, pxcrr, adj, hb , xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, rhoum,rhouc,rhouh, rhou,rhoud, rhouf, sigc, sigh, sigm, ss, iss_noa, rss_noa, yss_noa, kss_noa, css_noa,wss_noa,mss_noa,  hss_noa, ivss_noa, ggss_noa, lamdass_noa ;
  @#endif
  @#if roten==1 
  parameters tau, pxcrr, adj, hb , xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rot, rhoum,rhouc,rhouh, rhou,rhoud, rhouf, sigc, sigh, sigm, ss, iss_noa, rss_noa, yss_noa, kss_noa, css_noa,wss_noa,mss_noa,  hss_noa, ivss_noa, ggss_noa, lamdass_noa ;
  @#endif
 @#endif
 @#if pindex==1
  @#if roten==0 
  parameters index, tau, pxcrr, adj, hb , xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, rhoum,rhouc,rhouh, rhou,rhoud, rhouf, sigc, sigh, sigm, ss, iss_noa, rss_noa, yss_noa, kss_noa, css_noa,wss_noa,mss_noa,  hss_noa, ivss_noa, ggss_noa, lamdass_noa ;
  @#endif
 @#endif


@#endif

rho_rot_=rho_rot; 

@#if tech==1

//** Start of data_techzero==1 section **//
@#if data_techzero==1


@#if eshock==16
 
 rho_eta=aa0(j,17);
 xxx = aa0(j,18);
 xpai = aa0(j,19); 
 xy =aa0(j,20); 


 
 hb=aa0(j,21); % Internal habit formation param.
 
 
 sigc =aa0(j,23); 
 sigm = 2.56; 
 ss=aa0(j,24);
 phi = aa0(j,25);
 psi = aa0(j,26);


 rhou = aa0(j,27);
 rhoi = aa0(j,28);
 rhoud =0.5;     % Not used, irrelevant. 
 rhouh = aa0(j,29); 
 rhouc = aa0(j,30); 
 rhouf = aa0(j,31); 


 @#if ashock==1
 //alfa =0.1530682;
 alfa =aa0(j,32);
  @#endif

@#if ashock==0
alfa =alfa_n;
@#endif

@#if ashock==2
alfa =mode_adj*aa0(j,32);
@#endif

@#endif

pxcrr=pxcrr_;
adj=adj_;
rhoum =0.40;
gamma =0.99; % discount factor  
delta = 0.1000/4;  %capital depreciation rate
theta = 0.33;
notec=0;
epsilon = 0.1/4;
tau =0.18;
sigh=2.0;


@#if roten==0
 rho =aa0(j,22);  %Calve pricing parameter, rho=0 is flexible price
@#endif
@#if roten==1
 rot =((rho_rot_)*psi*(phi-1))/(  1-rho_rot_*psi-rho_rot_*(psi)*gamma + ((rho_rot_*psi)^2)*gamma ) ; % Rotember pricing parameter
@#endif



eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5)+(rho_eta^6)+(rho_eta^7))^(-1) );
eta2=eta1*(rho_eta);
eta3=eta1*(rho_eta)^2; 
eta4=eta1*(rho_eta)^3;
eta5=eta1*(rho_eta)^4;
eta6=eta1*(rho_eta)^5;
eta7=eta1*(rho_eta)^6;
eta8=eta1*(rho_eta)^7;


//** Steady state values: "ss" indicates staedy state.
iss = (1/gamma)-1;
rss = iss + delta;
wss = (((phi-1)/phi)*(1-theta)*(((1-theta)/theta)^(-theta))*(rss^(-theta)))^(1/(1-theta));
yss_hss = ((theta/(1-theta))*(wss/rss))^theta;
kss_hss = (yss_hss)^(1/theta);
nfss= (1-psi)/epsilon;


g1_t8 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8)  )^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
g2_t8 = (1-g1_t8)*(yss_hss) - (delta*kss_hss)- tau*(yss_hss);
g1 = g1_t8; 
g2 = g2_t8;
css_hss = g2_t8; 
etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8) );

hss= (wss*(css_hss^(-sigc))*((1-hb)^(-sigc))*(1-hb*gamma))^(1/(sigh+sigc)); 
css=(css_hss)*hss;
yss = (((theta/(1-theta))*(wss/rss))^theta)*hss;
kss = (yss_hss^(1/theta))*hss;

rdss = g1*yss;
x=((nfss)^(1+alfa))/rdss; % x is (kappa_bar*delta) in the paper

pfss = (1/(1-psi))*etas*rdss; % pfss is the steady state value of th epresent value of monopoly profits stream 
lamdass = (1-hb*gamma)*((1-hb)*css)^(-sigc); % lamdass = the steady state velu of CHI in the paper
mss = (((1-hb)*css)^(sigc/sigm))*(((1+iss)/iss)^(1/sigm))*((1-hb*gamma)^(-1/sigm));
ivss = delta*kss;
ggss =tau*yss; % ggss=staedy state value of gg, govt expnditure

@#endif
//** End of data_techzero==1 section **//


@#endif



@#if tech==0


@#if data_techzero==0  

@#if eshock==16


@#if pindex==0
 xxx = aa0(j,17);
 xpai = aa0(j,18);
 xy =aa0(j,19); 

 hb=aa0(j,20); % Internal habit formation param.
 sigc =aa0(j,22); 
 sigm = 2.56; 
 ss=aa0(j,23);
 phi = aa0(j,24);

 rhou = aa0(j,25);
 rhoi = aa0(j,26);
 rhoud = 0.5;
 rhouh = aa0(j,27);
 rhouc = aa0(j,28); 
 rhouf = aa0(j,29); 
@#endif

@#if pindex==1
 index=aa0(j,17); 
 xxx = aa0(j,18);
 xpai = aa0(j,19); 
 xy =aa0(j,20); 

 hb=aa0(j,21); % Internal habit formation param.
 sigc =aa0(j,23); 
 sigm = 2.56; 
 ss=aa0(j,24);
 phi = aa0(j,25);

 rhou = aa0(j,26);
 rhoi = aa0(j,27);
 rhoud = 0.5;
 rhouh = aa0(j,28); 
 rhouc = aa0(j,29); 
 rhouf = aa0(j,30); 
@#endif 




@#endif


pxcrr=pxcrr_;
adj=adj_;
rhoum =0.9;
sigh=2.0;
gamma =0.99; % discount factor  
delta = 0.10000000/4;  %capital depreciation rate
theta = 1/3;
notec=0;
tau=0.18;
//phi=6.0;

@#if pindex==0
 @#if roten==0
 rho =aa0(j,21);  //** rho=0 is flexible price
 @#endif
 @#if roten==1
 rot =((rho_rot_)*(phi-1))/(  1-rho_rot_-rho_rot_*gamma + ((rho_rot_)^2)*gamma ) ;
 @#endif
@#endif

@#if pindex==1
 @#if roten==0
 rho =aa0(j,22);  //** rho=0 is flexible price
 @#endif
@#endif


//** Steady state values: "ss" indicates staedy state.
iss_noa = (1/gamma)-1;
rss_noa = iss_noa + delta;
wss_noa = (((phi-1)/phi)*(1-theta)*(((1-theta)/theta)^(-theta))*(rss_noa^(-theta)))^(1/(1-theta));
yss_hss_noa = ((theta/(1-theta))*(wss_noa/rss_noa))^theta;
kss_hss_noa = (yss_hss_noa)^(1/theta);
css_hss_noa = yss_hss_noa - delta*kss_hss_noa- tau*yss_hss_noa; //* including the gov expenditure part "-tau*yss_hss_noa"
hss_noa = (wss_noa*(css_hss_noa^(-sigc))*((1-hb)^(-sigc))*(1-hb*gamma))^(1/(sigh+sigc));
yss_noa = (((theta/(1-theta))*(wss_noa/rss_noa))^theta)*hss_noa;
kss_noa = (yss_hss_noa^(1/theta))*hss_noa;
css_noa = css_hss_noa*hss_noa;
mss_noa = (((1-hb)*css_noa)^(sigc/sigm))*(((1+iss_noa)/iss_noa)^(1/sigm))*((1-hb*gamma)^(-1/sigm));
ivss_noa = delta*kss_noa;
ggss_noa =tau*yss_noa; % ggss=staedy state value of gg, govt expnditure
lamdass_noa = (1-hb*gamma)*((1-hb)*css_noa)^(-sigc);
@#endif

@#endif


//***************************************//
//*********** MODEL EQNS***************//
//************************************//
model(linear);

////////////////////////
//**** EPD-NK Model***//
///////////////////////

@#if tech==1

 @#if TI==1
rd=(1+alfa)*(eta1*nf(-7)+eta2*nf(-6)+eta3*nf(-5)+eta4*nf(-4)+eta5*nf(-3)+eta6*nf(-2)+eta7*nf(-1)+eta8*nf);

@#if roten==0
pf(+8)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*i(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*i(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*i(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*i(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*i(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*i(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*i(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*i + eta1*(1+iss)*(p(+7)-p(+8)) + eta2*((1+iss)^2)*(p(+6)-p(+8)) + eta3*((1+iss)^3)*(p(+5)-p(+8)) + eta4*((1+iss)^4)*(p(+4)-p(+8)) + eta5*((1+iss)^5)*(p(+3) -p(+8)) + eta6*((1+iss)^6)*(p(+2) -p(+8)) + eta7*((1+iss)^7)*(p(+1) -p(+8)) + eta8*((1+iss)^8)*(p -p(+8))  );    
@#endif
@#if roten==1
pf(+8)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*i(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*i(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*i(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*i(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*i(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*i(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*i(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*i + eta1*(1+iss)*(p(+7)-p(+8)) + eta2*((1+iss)^2)*(p(+6)-p(+8)) + eta3*((1+iss)^3)*(p(+5)-p(+8)) + eta4*((1+iss)^4)*(p(+4)-p(+8)) + eta5*((1+iss)^5)*(p(+3) -p(+8)) + eta6*((1+iss)^6)*(p(+2) -p(+8)) + eta7*((1+iss)^7)*(p(+1) -p(+8)) + eta8*((1+iss)^8)*(p -p(+8))  );    
@#endif



a = (1-psi)*nf(-7) + psi*a(-1);
 @#endif

 @#if TI==0
rd=(1+alfa)*(nf);
pf(+1)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* ( (i)  +(1+iss)*(p-p(+1)) );    
a = (1-psi)*nf + psi*a(-1);
 @#endif


c-hb*c(-1)= (1+gamma*hb)*(c(+1)-hb*c) - gamma*hb*(c(+2)-hb*c(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss/(1+iss))*(1/iss)*i + p - p(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*h = w - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(c - hb*c(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - uc(+1) - ud(+1) +ud );
sigm*(m-p) = (-1/(1+iss))*(1/iss)*i + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(c-hb*c(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) - ud) );
lamda = (1-delta)*gamma*lamda(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - (1/rss)*r(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(c(+2)-hb*c(+1)) - (1/rss)*r(+1) -uc(+2) -ud(+2) ) );
iv =  (gamma/(1+gamma))*iv(+1) + (1/(1+gamma))*iv(-1) + (1/(ss*(1+gamma)))*lamda  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(c-hb*c(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) ) ) ;
yss*y = css*c + ivss*iv + rdss*rd +ggss*gg;
(1-delta)*k(-1) +delta*iv - k = 0;
tiv= (ivss*iv + rdss*rd)/(ivss+rdss); //*total investment: R&D investment+ capital investment
y = (1/(phi-1))*a(-1) + u + theta*k(-1) + (1-theta)*h;
h - k(-1) - (1/rss)*r + w = 0;

@#if roten==0
pai = gamma*(pai(+1)) + (((1-gamma*rho*psi)*(1-rho*psi))/(rho*psi))*(-u + theta*(1/rss)*r +(1-theta)*w - (1/(phi-1))*a(-1)) + (gamma/(phi-1))*(a - a(-1)) - (1/(phi-1))*(a(-1) - a(-2));
@#endif
@#if roten==1
pai = gamma*psi*(pai(+1)) + ((phi-1)/rot)*(-u + theta*(1/rss)*r +(1-theta)*w ) + (gamma*psi/(phi-1))*(a) -(((phi-1)+rot+gamma*psi*rot)/((phi-1)*rot))* a(-1) + (1/(phi-1))*a(-2);
@#endif

@#if roten==0
pfss*pf - (gamma*psi)*pfss*pf(+1) + (gamma*psi)*pfss*(p - p(+1) + iss*(1/iss)*i) = (1/phi)*yss*y +((phi-1)/phi)*yss*u - ((phi-1)/phi)*(1-theta)*yss*w -  ((phi-1)/phi)*theta*yss*(1/rss)*r;
@#endif
@#if roten==1
pfss*pf - (gamma*psi)*pfss*pf(+1) = (1/phi)*yss*y +((phi-1)/phi)*yss*u - ((phi-1)/phi)*(1-theta)*yss*w -  ((phi-1)/phi)*theta*yss*(1/rss)*r;
@#endif

pai = p - p(-1);
gg=(1/tau)*uf +y; // * govt expenditure

y_g=y-y(-1);
c_g=c-c(-1);
tiv_g=tiv-tiv(-1);
h_g=h-h(-1);
w_g=w-w(-1);


/////////////////////////////////////////
// Technology and interest rate shocks //
////////////////////////////////////////

iss*(1/iss)*i = ((1+iss)^(xxx))*(xxx*(iss/(1+iss))*(1/iss)*i(-1) + (1-xxx)*((xpai)*(pai)+(xy)*(yg)) + v);

@#if ishock==1
 v = rhoi*v(-1) + ei;
@#endif

@#if ishock==0
 v = rhoi*v(-1);
@#endif

@#if fshock==1
uf = rhouf*uf(-1)+ef;
@#endif

@#if fshock==0
uf = rhouf*uf(-1);
@#endif


@#if ushock==1
 @#if news==1
    @#if eshock==16
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16); 
   @#endif 
   @#if eshock==20
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue16(-16)+eue17(-17)+eue18(-18)+eue19(-19)+eue(-20); 
   @#endif
 @#endif
 @#if news==0
   u = rhou*u(-1) +eu ; 
 @#endif
@#endif 

@#if ushock==0
 @#if news==1
  @#if eshock==16
   u = rhou*u(-1)  + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16); 
   @#endif 
   @#if eshock==20
   u = rhou*u(-1)  + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue16(-16)+eue17(-17)+eue18(-18)+eue19(-19)+eue(-20); 
   @#endif
 @#endif
 @#if news==0
 u = rhou*u(-1);
 @#endif 
@#endif

 
////////////////////////
// Preference shocks //
//////////////////////
@#if cshock==0
uc = rhouc*uc(-1); //* Preference shock
@#endif
@#if cshock==1
uc = rhouc*uc(-1) +euc; //* Preference shock
@#endif

@#if hshock==0
uh = rhouh*uh(-1); //* Preference shock
@#endif
@#if hshock==1
uh = rhouh*uh(-1) +euh; //* Preference shock
@#endif

@#if dshock==0
ud = rhoud*ud(-1);
@#endif
@#if dshock==1
ud = rhoud*ud(-1)+eud;
@#endif

um = rhoum*um(-1);



///////////////////////////
/////Others/////////////
//////////////////////

yg=y-ynobs; //* gdp gap 
hg = h-hn;

@#if roten==0
pinf= (gamma/(phi-1))*(a - a(-1)) - (1/(phi-1))*(a(-1) - a(-2))+(((1-gamma*rho*psi)*(1-rho*psi))/(rho*psi))*(-(1/(phi-1))*a(-1)) ;
@#endif
@#if roten==1
pinf=(gamma*psi/(phi-1))*(a) -(((phi-1)+rot+gamma*psi*rot)/((phi-1)*rot))* a(-1) + (1/(phi-1))*a(-1);
@#endif


mcgap=(-u + theta*(1/rss)*r +(1-theta)*w - (1/(phi-1))*an(-1)); //* (mc - mcn)
mcgape=mcgap+mee; 
mc=-u + theta*(1/rss)*r +(1-theta)*w; //*real marginal cost under sticky price
//mc=-(1/(phi-1))*a(-1)-u + theta*(1/rss)*r +(1-theta)*w; //*real marginal cost under sticky price
mce=mc+mee;
fmc=mc(+1);
pfmc=fmc(-1);
fpai = pai(+1);
ppai = pai(-1);
ppai2 = pai(-2);
ppai3 = pai(-3);
ppai4 = pai(-4);

fpaie = pxcrr*pai(+1)+meee;
pfpaie=fpaie(-1);
paiee=pai+(adj*gamma*meee)+(me);
paie=pai+me;//*pai with measument errors
ppaie=paie(-1);
shk_me=me;
shk_mee=mee;
shk_meee=meee;
ffp=p(+2);
fp=p(+1);
ffpai=ffp-fp;
pffpai=ffpai(-1);

@#if roten==0
inftr =  (((1-gamma*rho*psi)*(1-rho*psi))/(rho*psi))*(- (1/(phi-1))*a(-1)) + (gamma/(phi-1))*(a - a(-1)) - (1/(phi-1))*(a(-1) - a(-2)); //*the residural part of inflation eqn
@#endif
@#if roten==1
inftr =  (gamma*psi/(phi-1))*(a) -(((phi-1)+rot+gamma*psi*rot)/((phi-1)*rot))* a(-1) + (1/(phi-1))*a(-1); //*the residural part of inflation eqn
@#endif


/////TFP////
tfp=(1/(phi-1))*a(-1) + u ;
tfpu=u;


errmc=mee;

//////////////////////////
//* Naatural variables *//
////////////////////////
  @#if TI==1
 rdn=(1+alfa)*(eta1*nfn(-7)+eta2*nfn(-6)+eta3*nfn(-5)+eta4*nfn(-4)+eta5*nfn(-3)+eta6*nfn(-2)+eta7*nfn(-1)+eta8*nfn);
 pfn(+8)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*in(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*in(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*in(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*in(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*in(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*in(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*in(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*in + eta1*(1+iss)*(pn(+7)-pn(+8)) + eta2*((1+iss)^2)*(pn(+6)-pn(+8)) + eta3*((1+iss)^3)*(pn(+5)-pn(+8)) + eta4*((1+iss)^4)*(pn(+4)-pn(+8)) + eta5*((1+iss)^5)*(pn(+3) -pn(+8)) + eta6*((1+iss)^6)*(pn(+2) -pn(+8)) + eta7*((1+iss)^7)*(pn(+1) -pn(+8)) + eta8*((1+iss)^8)*(pn -pn(+8))  );    
 an = (1-psi)*nfn(-7) + psi*an(-1);
   @#endif

 @#if TI==0
rdn=(1+alfa)*(nfn);
pfn(+1)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* ( (in)  +(1+iss)*(pn-pn(+1)) );    
an = (1-psi)*nfn + psi*an(-1);
  @#endif

cn-hb*cn(-1)= (1+gamma*hb)*(cn(+1)-hb*cn) - gamma*hb*(cn(+2)-hb*cn(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss/(1+iss))*(1/iss)*in + pn - pn(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*hn = wn - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(cn - hb*cn(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - uc(+1) - ud(+1) +ud );
sigm*(0-pn) = (-1/(1+iss))*(1/iss)*in + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(cn-hb*cn(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) - ud) );
lamdan = (1-delta)*gamma*lamdan(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - (1/rss)*rn(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(cn(+2)-hb*cn(+1)) - (1/rss)*rn(+1) -uc(+2) -ud(+2) ) );
ivn =  (gamma/(1+gamma))*ivn(+1) + (1/(1+gamma))*ivn(-1) + (1/(ss*(1+gamma)))*lamdan  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(cn-hb*cn(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) ) ) ;
//yss*yn = css*cn + ivss*ivn + rdss*rdn;
yss*yn = css*cn + ivss*ivn + rdss*rdn +ggss*ggn;
(1-delta)*kn(-1) +delta*ivn - kn = 0;
tivn= (ivss*ivn + rdss*rdn)/(ivss+rdss);
ynobs=yn;
hn - kn(-1) - (1/rss)*rn + wn = 0;
yn = (1/(phi-1))*an(-1) + u + theta*kn(-1) + (1-theta)*hn;
0 = (((1-gamma*0*psi)*(1-0*psi)))*(-u + theta*(1/rss)*rn +(1-theta)*wn - (1/(phi-1))*an(-1)) ;

@#if roten==0
pfss*pfn - (gamma*psi)*pfss*pfn(+1) + (gamma*psi)*pfss*(pn - pn(+1) + iss*(1/iss)*in) = (1/phi)*yss*yn +((phi-1)/phi)*yss*u - ((phi-1)/phi)*(1-theta)*yss*wn -  ((phi-1)/phi)*theta*yss*(1/rss)*rn;
@#endif
@#if roten==1
pfss*pfn - (gamma*psi)*pfss*pfn(+1)  = (1/phi)*yss*yn +((phi-1)/phi)*yss*u - ((phi-1)/phi)*(1-theta)*yss*wn -  ((phi-1)/phi)*theta*yss*(1/rss)*rn;
@#endif


pain = pn - pn(-1);

tfpn=(1/(phi-1))*an(-1) + u;

ggn= (1/tau)*uf +yn;

end;

initval;
ynobs=0;
tiv=0;
p = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;

a = 0;
i=0;
lamda =0;
iv=0;

v=0;


nf=0;
nfn=0;


y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;
rd = 0;
pf =0;
//q=0;
//eue=0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
//mn = 0;
wn =0 ;
rn = 0;
rdn = 0;
pfn =0;
pn = 0; 
pain = 0;
an = 0;
in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;

end;

steady_state_model;
ynobs=0;
tiv=0;
p = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;

a = 0;
i=0;
lamda =0;


v=0;


nf=0;
nfn=0;


y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;
rd = 0;
pf =0;
//q=0;
//eue=0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
//mn = 0;
wn =0 ;
rn = 0;
rdn = 0;
pfn =0;
pn = 0; 
pain = 0;
an = 0;
in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;

end;

@#endif


/////////////////////////////////////////
//**** Basic and Standard NK models***//
///////////////////////////////////////

@#if tech==0

c-hb*c(-1)= (1+gamma*hb)*(c(+1)-hb*c) - gamma*hb*(c(+2)-hb*c(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss_noa/(1+iss_noa))*(1/iss_noa)*i + p - p(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*h = w - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(c - hb*c(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - uc(+1) - ud(+1) +ud );
sigm*(m-p) = (-1/(1+iss_noa))*(1/iss_noa)*i + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(c-hb*c(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) - ud) );
lamda = (1-delta)*gamma*lamda(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - (1/rss_noa)*r(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(c(+2)-hb*c(+1)) - (1/rss_noa)*r(+1) -uc(+2) -ud(+2) ) );
iv =  (gamma/(1+gamma))*iv(+1) + (1/(1+gamma))*iv(-1) + (1/(ss*(1+gamma)))*lamda  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(c-hb*c(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) ) ) ;


yss_noa*y = css_noa*c + ivss_noa*iv ;
//** even with rotember pricing this does not change **//


(1-delta)*k(-1) +delta*iv - k = 0;

y = u + theta*k(-1) + (1-theta)*h;
h - k(-1) - (1/rss_noa)*r + w = 0;

//pai = gamma*(pai(+1)) + (((1-gamma*rho)*(1-rho))/(rho))*(-u + theta*(1/rss_noa)*r +(1-theta)*w ) ;
@#if pindex==0
 
 @#if roten==0
 pai = gamma*(pai(+1)) + (((1-gamma*rho)*(1-rho))/(rho))*(-u + theta*(1/rss_noa)*r +(1-theta)*w ) ;
 @#endif

 @#if roten==1
 pai = gamma*(pai(+1)) + ((phi-1)/rot)*(-u + theta*(1/rss_noa)*r +(1-theta)*w ) ;
 @#endif

@#endif

@#if pindex==1

 @#if roten==0
 pai = ( gamma/(1+gamma*rho*index) )*(pai(+1)) + ( index/(1+gamma*rho*index) )*pai(-1)+ (  ((1-gamma*rho)*(1-rho))/( (rho)*(1+gamma*rho*index) )  )*(-u + theta*(1/rss_noa)*r +(1-theta)*w ); 
 @#endif
 

@#endif

pai = p - p(-1);

tq = sigc*c+lamda; //*Tobin's Q
//tq = -(css_noa^(-sigc))*sigc*c-(1/lamdass_noa)*lamda; //*Tobin's Q

gg=(1/tau)*uf +y; // * govt expenditure

y_g=y-y(-1);
c_g=c-c(-1);
iv_g=iv-iv(-1);
h_g=h-h(-1);
w_g=w-w(-1);


//////////////////////////////////////////
// Technology and interest rate shocks //
////////////////////////////////////////
iss_noa*(1/iss_noa)*i = ((1+iss_noa)^(xxx))*(xxx*(iss_noa/(1+iss_noa))*(1/iss_noa)*i(-1) + (1-xxx)*((xpai)*(pai)+(xy)*(yg)) + v);

@#if ishock==1
 v = rhoi*v(-1) + ei;
@#endif

@#if ishock==0
 v = rhoi*v(-1);
@#endif

@#if fshock==1
uf = rhouf*uf(-1)+ef;
@#endif

@#if fshock==0
uf = rhouf*uf(-1);
@#endif

@#if ushock==1
 @#if news==1
   @#if eshock==16
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16); 
   @#endif 
   @#if eshock==20
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue16(-16)+eue17(-17)+eue18(-18)+eue19(-19)+eue(-20); 
   @#endif
 @#endif
 @#if news==0
    u = rhou*u(-1);
 @#endif
@#endif 

@#if ushock==0
 @#if news==1
   @#if eshock==16
   u = rhou*u(-1) + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16); 
   @#endif 
   @#if eshock==20
   u = rhou*u(-1)  + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue16(-16)+eue17(-17)+eue18(-18)+eue19(-19)+eue(-20); 
   @#endif
 @#endif

 @#if news==0
 u = rhou*u(-1);
 @#endif 
@#endif

 
////////////////////////
// Preference shocks //
//////////////////////
@#if dshock==1
ud = rhoud*ud(-1) +eud; //* Preference shock
@#endif

@#if dshock==0
ud = rhoud*ud(-1); //* Preference shock
@#endif

@#if cshock==1
uc = rhouc*uc(-1) +euc; //* Preference shock
@#endif

@#if cshock==0
uc = rhouc*uc(-1); //* Preference shock
@#endif


@#if hshock==1
uh = rhouh*uh(-1) +euh; //* Preference shock
@#endif

@#if hshock==0
uh = rhouh*uh(-1); //* Preference shock
@#endif


//um = rhoum*um(-1) +eum; //* Preference shock
um = rhoum*um(-1) ; //* Preference shock

///////////////////////////
/////Others/////////////
//////////////////////


yg=y-ynobs; //* gdp gap 
hg = h-hn;


mcgap=(-u + theta*(1/rss_noa)*r +(1-theta)*w ); //* deviation of marginal costs from its flwxible level
mcgape=mcgap+mee; //* maginal cost deviationa with measuremant errors
mc=mcgap; //* deviation of marginal costs from its flwxible level
mce=mcgape; //* maginal cost deviationa with measuremant errors
fmc=mc(+1);
pfmc=fmc(-1);
fpai = pai(+1);
ppai = pai(-1);

ppai2 = pai(-2);
ppai3 = pai(-3);
ppai4 = pai(-4);

fpaie = pxcrr*pai(+1)+meee;
pfpaie=fpaie(-1);
paiee=pai+(adj*gamma*meee)+(me);//*pai with 2sls errors & measument errors
paie=pai+me;//*pai with measument errors
ppaie=paie(-1);
shk_me=me;
shk_mee=mee;
shk_meee=meee;

ffp=p(+2);
fp=p(+1);
ffpai=ffp-fp;
pffpai=ffpai(-1);

tfp=u;
tfpu=u;

errmc=mee;
//////////////////////////
//* natural variables *//
////////////////////////


cn-hb*cn(-1)= (1+gamma*hb)*(cn(+1)-hb*cn) - gamma*hb*(cn(+2)-hb*cn(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss_noa/(1+iss_noa))*(1/iss_noa)*in + pn - pn(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*hn = wn - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(cn - hb*cn(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - uc(+1) - ud(+1) +ud );
sigm*(0-pn) = (-1/(1+iss_noa))*(1/iss_noa)*in + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(cn-hb*cn(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) - ud) );

lamdan = (1-delta)*gamma*lamdan(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - (1/rss_noa)*rn(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(cn(+2)-hb*cn(+1)) - (1/rss_noa)*rn(+1) -uc(+2) -ud(+2) ) );
ivn =  (gamma/(1+gamma))*ivn(+1) + (1/(1+gamma))*ivn(-1) + (1/(ss*(1+gamma)))*lamdan  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(cn-hb*cn(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) ) ) ;

yss_noa*yn = css_noa*cn + ivss_noa*ivn +ggss_noa*ggn;

(1-delta)*kn(-1) +delta*ivn - kn = 0;
ynobs=yn;

hn - kn(-1) - (1/rss_noa)*rn + wn = 0;
yn = u + theta*kn(-1) + (1-theta)*hn;
0 = (((1-gamma*0)*(1-0)))*(-u + theta*(1/rss_noa)*rn +(1-theta)*wn ) ;

pain = pn - pn(-1);

tfpn=u;

ggn= (1/tau)*uf +yn;
end;


initval;

p = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;


i=0;
lamda =0;
iv=0;


v=0;



y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
wn =0 ;
rn = 0;
//rdn = 0;

pn = 0; 
pain = 0;

in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;

end;

steady_state_model;
p = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;


i=0;
lamda =0;
//tiv=0;

v=0;




y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;


yn =0 ;
kn = 0;
hn = 0;
cn = 0;
wn =0 ;
rn = 0;

pn = 0; 
pain = 0;

in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;
end;
@#endif


 shocks;


 var me;  // Not used: measurement errors  for 'mc'
 stderr me_std; 

 var mee;  // Not Used: measurement errors  for 'mc'
 stderr mee_std;

 var meee;  // Not used 
 stderr meee_std;

//***********************************//
//**** Start of fig==0 Section **////
//***********************************//
@#if fig==0
@#if data_techzero==0  
 
 var eu;
 stderr aa0(j,1); 
 
 var ei;
 stderr aa0(j,2);

   var euc;
 stderr aa0(j,3);

  var euh;
 stderr aa0(j,4);

 var ef;
 stderr aa0(j,5);

@#endif

@#if data_techzero==1  
 
   var eu;
 stderr aa0(j,1); 
 
 var ei;
 stderr aa0(j,2);

  var euc;
 stderr aa0(j,3);

  var euh;
 stderr aa0(j,4);

  var ef;
 stderr aa0(j,5);

@#endif
   

///////////////////
// NEWS SHOCKS//
//////////////////
@#if data_techzero==0 
 
var eue1;
 stderr 0;

  var eue2;
 stderr 0;

  var eue3;
 stderr 0;

  var eue4;
 stderr aa0(j,6);

  var eue5;
 stderr 0;

  var eue6;
 stderr 0;

  var eue7;
 stderr  0;

 var eue8;
 stderr aa0(j,7);

  var eue9;
 stderr  0;
 
  var eue10;
 stderr  0;
   
 @#if eshock==16
 
 var eue11;
 stderr 0;

  var eue12;
 stderr aa0(j,8);

 var eue13;
 stderr 0;
 
 var eue14;
 stderr 0;

 var eue15;
 stderr 0;
 
 var eue;
 stderr aa0(j,9);


corr eue4, eue8 =aa0(j,11); corr eue4, eue12 =aa0(j,12); corr eue4, eue =aa0(j,13); corr eue4, eu =0;
corr eue8, eue12= aa0(j,14); corr eue8, eue =aa0(j,15); corr eue8, eu =0;
corr eue12, eue =aa0(j,16); corr eue12, eu =0;
corr eue, eu =0;

@#endif


@#endif

@#if data_techzero==1 

var eue1;
 stderr 0;

  var eue2;
 stderr 0;

  var eue3;
 stderr 0;

  var eue4;
 stderr aa0(j,6);

  var eue5;
 stderr 0;

  var eue6;
 stderr 0;

  var eue7;
 stderr  0;

 var eue8;
 stderr aa0(j,7);

  var eue9;
 stderr  0;
 
  var eue10;
 stderr  0;
  
 
 @#if eshock==16
 var eue11;
 stderr  0;

  var eue12;
 stderr aa0(j,8);

 var eue13;
 stderr  0;
 
 var eue14;
 stderr  0;

 var eue15;
 stderr  0;
 
 var eue;
 stderr aa0(j,9);

corr eue4, eue8 =aa0(j,11); corr eue4, eue12 =aa0(j,12); corr eue4, eue =aa0(j,13); corr eue4, eu =0;
corr eue8, eue12 = aa0(j,14); corr eue8, eue =aa0(j,15); corr eue8, eu =0;
corr eue12, eue =aa0(j,16); corr eue12, eu =0;
corr eue, eu =0;
 
  @#endif


@#endif
@#endif
//***********************************//
//** End of fig==0 section **//
//***********************************//


//***********************************//
//** Start fig==1 section **//
//***********************************//
@#if fig==1

@#if cndvar==0

@#if data_techzero==0  
 
 @#if shki==1
 var ei;
 stderr 0.01;
 @#endif
 @#if shki==0
 var ei;
 stderr 0;
 @#endif

 @#if shku==1
 var eu;
 stderr 0.01; 
 @#endif
 @#if shku==0
 var eu;
 stderr 0; 
 @#endif

 
 @#if shkuc==1
 var euc;
 stderr 0.01;
 @#endif
 @#if shkuc==0
 var euc;
 stderr 0;
 @#endif

 @#if shkuh==1
 var euh;
 stderr 0.01;
 @#endif
 @#if shkuh==0
 var euh;
 stderr 0;
 @#endif

@#if shkf==1
 var ef;
 stderr 0.01;
 @#endif
 @#if shkf==0
 var ef;
 stderr 0;
 @#endif 
 
@#endif

@#if data_techzero==1  
 
 @#if shki==1
 var ei;
 stderr 0.01;
 @#endif
 @#if shki==0
 var ei;
 stderr 0;
 @#endif
 
 @#if shku==1
 var eu;
 stderr 0.01; 
 @#endif
 @#if shku==0
 var eu;
 stderr 0; 
 @#endif

  @#if shkuc==1
 var euc;
 stderr 0.01;
 @#endif
 @#if shkuc==0
 var euc;
 stderr 0;
 @#endif

 @#if shkuh==1
 var euh;
 stderr 0.01;
 @#endif
 @#if shkuh==0
 var euh;
 stderr 0;
 @#endif 

@#if shkf==1
 var ef;
 stderr 0.01;
 @#endif
 @#if shkf==0
 var ef;
 stderr 0;
 @#endif 

@#endif



///////////////////
// NEWS SHOCKS//
//////////////////
 @#if shkue1==1
  var eue1;
 stderr 0.01;
 @#endif
  @#if shkue1==0
  var eue1=0;
 @#endif
 
 @#if shkue2==1
  var eue2;
 stderr 0.01;
 @#endif
 @#if shkue2==0
  var eue2=0;
 @#endif
 
 @#if shkue3==1
  var eue3;
 stderr 0.01;
 @#endif
 @#if shkue3==0
  var eue3=0;
 @#endif

 @#if shkue4==1
  var eue4;
 stderr 0.01;
 @#endif
  @#if shkue4==0
  var eue4=0;
 @#endif

@#if shkue5==1
  var eue5;
 stderr 0.01;
 @#endif
 
@#if shkue5==0
  var eue5=0;
 @#endif

@#if shkue6==1
  var eue6;
 stderr 0.01;
 @#endif
 
 @#if shkue6==0
  var eue6=0;
 @#endif

@#if shkue7==1
  var eue7;
 stderr 0.01;
 @#endif
 
 @#if shkue7==0
  var eue7=0;
 @#endif

@#if shkue8==1
  var eue8;
 stderr 0.01;
 @#endif
 
 @#if shkue8==0
  var eue8=0;
 @#endif

@#if shkue9==1
  var eue9;
 stderr 0.01;
 @#endif
 
 @#if shkue9==0
  var eue9=0;
 @#endif

 @#if eshock==10
  @#if shkue==1
  var eue;
  stderr 0.01;
  @#endif
 
  @#if shkue==0
  var eue=0;
  @#endif
 @#endif

 @#if eshock==12
  @#if shkue10==1
  var eue10;
  stderr 0.01;
  @#endif
 
  @#if shkue10==0
  var eue10=0;
  @#endif
  
  @#if shkue11==1
  var eue11;
  stderr 0.01;
  @#endif
 
  @#if shkue11==0
  var eue11=0;
  @#endif
   
  @#if shkue==1
  var eue;
  stderr 0.01;
  @#endif
 
  @#if shkue==0
  var eue=0;
  @#endif
 @#endif



@#if eshock==16
  @#if shkue10==1
  var eue10;
  stderr 0.01;
  @#endif
 
  @#if shkue10==0
  var eue10=0;
  @#endif
  
  @#if shkue11==1
  var eue11;
  stderr 0.01;
  @#endif
 
  @#if shkue11==0
  var eue11=0;
  @#endif
   
  @#if shkue12==1
  var eue12;
  stderr 0.01;
  @#endif
 
  @#if shkue12==0
  var eue12=0;
  @#endif

  @#if shkue13==1
  var eue13;
  stderr 0.01;
  @#endif
 
  @#if shkue13==0
  var eue13=0;
  @#endif
  

  @#if shkue14==1
  var eue14;
  stderr 0.01;
  @#endif
 
  @#if shkue14==0
  var eue14=0;
  @#endif

  @#if shkue15==1
  var eue15;
  stderr 0.01;
  @#endif
 
  @#if shkue15==0
  var eue15=0;
  @#endif


   @#if shkue==1
  var eue;
  stderr 0.01;
  @#endif
 
  @#if shkue==0
  var eue=0;
  @#endif
 @#endif

//@#endif




@#if eshock==20
  @#if shkue10==1
  var eue10;
  stderr 0.01;
  @#endif
 
  @#if shkue10==0
  var eue10=0;
  @#endif
  
  @#if shkue11==1
  var eue11;
  stderr 0.01;
  @#endif
 
  @#if shkue11==0
  var eue11=0;
  @#endif
   
  @#if shkue12==1
  var eue12;
  stderr 0.01;
  @#endif
 
  @#if shkue12==0
  var eue12=0;
  @#endif

  @#if shkue13==1
  var eue13;
  stderr 0.01;
  @#endif
 
  @#if shkue13==0
  var eue13=0;
  @#endif
  

  @#if shkue14==1
  var eue14;
  stderr 0.01;
  @#endif
 
  @#if shkue14==0
  var eue14=0;
  @#endif

  @#if shkue15==1
  var eue15;
  stderr 0.01;
  @#endif
 
  @#if shkue15==0
  var eue15=0;
  @#endif


   @#if shkue16==1
  var eue16;
  stderr 0.01;
  @#endif
 
  @#if shkue16==0
  var eue16=0;
  @#endif

  @#if shkue17==1
  var eue17;
  stderr 0.01;
  @#endif
 
  @#if shkue17==0
  var eue17=0;
  @#endif
  

  @#if shkue18==1
  var eue18;
  stderr 0.01;
  @#endif
 
  @#if shkue18==0
  var eue18=0;
  @#endif

  @#if shkue19==1
  var eue19;
  stderr 0.01;
  @#endif
 
  @#if shkue19==0
  var eue19=0;
  @#endif

   @#if shkue==1
  var eue;
  stderr 0.01;
  @#endif
 
  @#if shkue==0
  var eue=0;
  @#endif
 @#endif

@#endif




@#if std_spec==1

@#if data_techzero==0  
 @#if shku==1
 var eu;
 stderr aa0(j,1); 
 @#endif
 @#if shku==0
 var eu;
 stderr 0; 
 @#endif

 @#if shki==1
 var ei;
 stderr aa0(j,2);
 @#endif
 @#if shki==0
 var ei;
 stderr 0;
 @#endif

 @#if shkuc==1
 var euc;
 stderr aa0(j,3);
 @#endif
 @#if shkuc==0
 var euc;
 stderr 0;
 @#endif

 @#if shkuh==1
 var euh;
 stderr aa0(j,4);
 @#endif
 @#if shkuh==0
 var euh;
 stderr 0;
 @#endif 

@#if shkf==1
 var ef;
 stderr aa0(j,5);
 @#endif
 @#if shkf==0
 var ef;
 stderr 0;
 @#endif 

@#endif

@#if data_techzero==1  
 @#if shku==1
 var eu;
 stderr aa0(j,1); 
 @#endif
 @#if shku==0
 var eu;
 stderr 0; 
 @#endif

 @#if shki==1
 var ei;
 stderr aa0(j,2);
 @#endif
 @#if shki==0
 var ei;
 stderr 0;
 @#endif

 @#if shkuc==1
 var euc;
 stderr aa0(j,3);
 @#endif
 @#if shkuc==0
 var euc;
 stderr 0;
 @#endif

 @#if shkuh==1
 var euh;
 stderr aa0(j,4);
 @#endif
 @#if shkuh==0
 var euh;
 stderr 0;
 @#endif 

@#if shkf==1
 var ef;
 stderr aa0(j,5);
 @#endif
 @#if shkf==0
 var ef;
 stderr 0;
 @#endif 

@#endif


///////////////////
// NEWS SHOCKS//
//////////////////
@#if data_techzero==0 
 
var eue1;
 stderr aa0(j,6);

  var eue2;
 stderr aa0(j,7);

  var eue3;
 stderr aa0(j,8);

  var eue4;
 stderr aa0(j,9);

  var eue5;
 stderr aa0(j,10);

  var eue6;
 stderr aa0(j,11);

  var eue7;
 stderr aa0(j,12);

 var eue8;
 stderr aa0(j,13);

  var eue9;
 stderr aa0(j,14);


@#if eshock==16
  
  var eue10;
 stderr aa0(j,15);
  
 var eue11;
 stderr aa0(j,16);

  var eue12;
 stderr aa0(j,17);

 var eue13;
 stderr aa0(j,18);
 
 var eue14;
 stderr aa0(j,19);

 var eue15;
 stderr aa0(j,20);
 
 var eue;
 stderr aa0(j,21);

corr eue4, eue8 =aa0(j,35); corr eue4, eue12 =aa0(j,36); corr eue4, eue =aa0(j,37); corr eue4, eu =aa0(j,38);
corr eue8, eue12= aa0(j,39); corr eue8, eue =aa0(j,40); corr eue8, eu =aa0(j,41);
corr eue12, eue =aa0(j,42); corr eue12, eu =aa0(j,43);
corr eue, eu =aa0(j,44);
 
 @#endif
 

@#if eshock==20
  var eue10;
 stderr aa0(j,14);
  
 var eue11;
 stderr aa0(j,15);

  var eue12;
 stderr aa0(j,16);

 var eue13;
 stderr aa0(j,17);

  var eue14;
 stderr aa0(j,18);

 var eue15;
 stderr aa0(j,19);
 
 var eue16;
 stderr aa0(j,20);

 var eue17;
 stderr aa0(j,21);

  var eue18;
 stderr aa0(j,22);

 var eue19;
 stderr aa0(j,23);
 
 var eue;
 stderr aa0(j,24);
 
 @#endif
 



@#endif

@#if data_techzero==1 
var eue1;
 stderr aa0(j,6);

  var eue2;
 stderr aa0(j,7);

  var eue3;
 stderr aa0(j,8);

  var eue4;
 stderr aa0(j,9);

  var eue5;
 stderr aa0(j,10);

  var eue6;
 stderr aa0(j,11);

  var eue7;
 stderr aa0(j,12);

 var eue8;
 stderr aa0(j,13);

  var eue9;
 stderr aa0(j,14);


@#if eshock==16
 var eue10;
 stderr aa0(j,15);
  
 var eue11;
 stderr aa0(j,16);

  var eue12;
 stderr aa0(j,17);

 var eue13;
 stderr aa0(j,18);
 
 var eue14;
 stderr aa0(j,19);

 var eue15;
 stderr aa0(j,20);
 
 var eue;
 stderr aa0(j,21);

corr eue4, eue8 =aa0(j,38); corr eue4, eue12 =aa0(j,39); corr eue4, eue =aa0(j,40); corr eue4, eu =aa0(j,41);
corr eue8, eue12 = aa0(j,42); corr eue8, eue =aa0(j,43); corr eue8, eu =aa0(j,44);
corr eue12, eue =aa0(j,45); corr eue12, eu =aa0(j,46);
corr eue, eu =aa0(j,47);


 @#endif
 

@#if eshock==20
   var eue10;
 stderr aa0(j,14);
  
 var eue11;
 stderr aa0(j,15);

  var eue12;
 stderr aa0(j,16);

 var eue13;
 stderr aa0(j,17);

  var eue14;
 stderr aa0(j,18);

 var eue15;
 stderr aa0(j,19);
 
 var eue16;
 stderr aa0(j,20);

 var eue17;
 stderr aa0(j,21);

  var eue18;
 stderr aa0(j,22);

 var eue19;
 stderr aa0(j,23);
 
 var eue;
 stderr aa0(j,24);
 
 @#endif
 

@#endif

@#endif


@#endif
//***********************************//
//** End of fig==1 section **//
//***********************************//
 

 
end;



for ii=iia:iib
set_dynare_seed(ii); 

@#if large_sample==1

@#if tech==1
 @#if fig==0
 stoch_simul(irf=30, order=1,periods=15000, nograph,noprint,qz_criterium=0.999999)  y yn a an c cn iv ivn k kn h hn tiv rd pai w r i rdn tiv mc mcgap nf nfn inftr;    
 @#endif

 @#if fig==1
 @#if cndvar==0
 stoch_simul(irf=30, order=1,periods=15000,nograph) y yn a an c cn iv ivn k kn h hn tiv rd pai w r i rdn tiv mc mcgap nf nfn inftr;  
 @#endif 
 @#if cndvar==1
 stoch_simul(irf=30, order=1,periods=15000,nograph, drop=300,conditional_variance_decomposition = [1 2 4 8 16 40]) y yn a an c cn iv ivn k kn h hn tiv rd pai w r i rdn tiv mc mcgap nf nfn inftr;  
 @#endif
@#endif
@#endif

@#if tech==0
 @#if fig==0
 stoch_simul(irf=30, order=1,periods=15000,nograph,noprint ,qz_criterium=0.999999)  y yn c cn iv ivn k kn h hn  pai w r i mc mcgap ; 
 @#endif   
 @#if fig==1
 @#if cndvar==0
 stoch_simul(irf=30, order=1,periods=15000,nograph) y yn c cn iv ivn k kn h hn  pai w r i mc mcgap ; 
 @#endif
 @#if cndvar==1
 stoch_simul(irf=30, order=1,periods=0,nograph, drop=300,conditional_variance_decomposition = [1 2 4 8 16 40]) y yn c cn iv ivn k kn h hn  pai w r i mc mcgap ; 
 @#endif
 @#endif   
@#endif

@#endif


@#if large_sample==0

@#if tech==1
 @#if fig==0
 stoch_simul(irf=30, order=1,periods=2000,nograph, noprint ,qz_criterium=0.999999)  y yn a an c cn iv ivn k kn h hn tiv rd pai w r i rdn tivn mc mcgap nf nfn inftr;    
 @#endif
 @#if fig==1
 @#if cndvar==0
 stoch_simul(irf=30, order=1,periods=2000,nograph ,qz_criterium=0.999999) y yn a an c cn iv ivn k kn h hn tiv rd pai w r i rdn tivn mc mcgap nf nfn inftr;  
 @#endif 
 @#if cndvar==1
 stoch_simul(irf=30, order=1,periods=0,nograph, drop=10,conditional_variance_decomposition = [1 2 4 8 16 40]) y yn a an c cn iv ivn k kn h hn tiv rd pai w r i rdn tivn mc mcgap nf nfn inftr;  
 @#endif
@#endif
@#endif

@#if tech==0
 @#if fig==0
 stoch_simul(irf=30, order=1,periods=2000,nograph,noprint ,qz_criterium=0.999999)  y yn c cn iv ivn k kn h hn  pai w r i mc mcgap ; 
 @#endif   
 @#if fig==1
 @#if cndvar==0
 stoch_simul(irf=30, order=1,periods=2000,nograph ,qz_criterium=0.999999) y yn c cn iv ivn k kn h hn  pai w r i mc mcgap ; 
 @#endif
 @#if cndvar==1
 stoch_simul(irf=30, order=1,periods=0,nograph, drop=300,conditional_variance_decomposition = [1 2 4 8 16 40]) y yn c cn iv ivn k kn h hn  pai w r i mc mcgap ; 
 @#endif
 @#endif   
@#endif

@#endif



//** start of fig==0 **//
@#if fig==0
//********************//


//sim_period=160;

cut_period1=cut_period+1;
cut_period2=cut_period+2;
cut_period3=cut_period+3;

//*****************************************************************************************************************************//
//** PC curve estimates: inflation = coef2*fpai(expected inflation) + coef3*MC(maginal costs) + coef4*ppai(lagged inflation) **//
//*******************************************************************************************************************************//
sim_period_=cut_period+sim_period-2;
pai_=pai(cut_period2:sim_period_+1);
paipls_=pai(cut_period3:sim_period_+2);
paie_=paie(cut_period3:sim_period_+1);
paiee_=paiee(cut_period2:sim_period_+2);

@#if scl_adj==0

 @#if corr_spec_mc==0
mce_=mce(cut_period1:sim_period_);
@#endif
@#if corr_spec_fpai==0
fpaie_=fpaie(cut_period1:sim_period_);
@#endif


@#if corr_spec_mc==1
//mce_=randwithcorr(mc(cut_period2:sim_period_+1), crr_mc, 0, stde_mc); 
//mce2_=randwithcorr(mc(cut_period2:sim_period_+1), mc_adj, 0, stde_mc); 
//fmce_=randwithcorr(pfmc(cut_period2:sim_period_+1), mc_adj2, 0, stde_fmc); 

mce_=randwithcorr(mc(cut_period2:sim_period_+1), crr_mc, 0, std(mc(cut_period2:sim_period_+1))); 
mce2_=randwithcorr(mc(cut_period2:sim_period_+1), mc_adj, 0, std(mc(cut_period2:sim_period_+1))); 
fmce_=randwithcorr(pfmc(cut_period2:sim_period_+1), mc_adj2, 0, std(pfmc(cut_period2:sim_period_+1))); 


@#endif
@#if corr_spec_fpai==1
//fpaie_=randwithcorr(fpai(1:sim_period_), crr_fpai, 0, stde_fpai);
//ffpaie_=randwithcorr(pffpai(cut_period2:sim_period_+1), crr_fpai, 0, stde_fpai);
ffpaie_=randwithcorr(pffpai(cut_period2:sim_period_+1), crr_fpai, 0, std(pffpai(cut_period2:sim_period_+1)));

@#endif

@#endif

@#if scl_adj==1
mce_=scl1*mce(cut_period2:sim_period_+1);
fpaie_=scl2*fpaie(cut_period2:sim_period_)+1;
@#endif

ppai_=ppai(cut_period2:sim_period_+1);
ppai2_=ppai2(cut_period2:sim_period_+1);
ppai3_=ppai3(cut_period2:sim_period_+1);
ppai4_=ppai4(cut_period2:sim_period_+1);

//X=[fpai(1:sim_period_), mcgap(1:sim_period_), ppai(1:sim_period_)];
//X=[fpai(1:sim_period_), mcgape(1:sim_period_), ppai(1:sim_period_)];

//tpaiepls_=fpaie((1:sim_period_)).*paipls_;
//tpai_=fpaie(1:sim_period_).*pai_;
//tmce_=fpaie(1:sim_period_).*mce(1:sim_period_);
//tppai_=fpaie(1:sim_period_).* ppaie(1:sim_period_);


@#if merror==7
sim_period__=sim_period_-cut_period;
@#if just_id==1
//***********IV estimation(just identified)*******************//
@#if nocon==0
x1_ = ones(sim_period__,1);
X_= [x1_ paipls_ mce_ ppai_];
 @#if mc_adj==0
 Z_= [x1_ fpaie_ mce_ ppai_];
 @#endif
 @#if mc_adj==1
 Z_= [x1_ fpaie_ mce2_ ppai_];
 @#endif
 @#if mc_adj==2
 Z_= [x1_ fpaie_ fmce_ ppai_];
 @#endif
@#endif

@#if nocon==1
X_= [ paipls_ mce_ ppai_];
Z_= [ fpaie_ mce_ ppai_];
@#endif
Y_ = pai_;
beta_hat_iv = (Z_'*X_)^(-1)*Z_'*Y_;    % IV estimated coeff
yhat = X_*beta_hat_iv;
res = Y_ - yhat; %residual ;
resvar = (res'*res)./(sim_period__);  %Var(residuals)
varcoef = resvar * inv(Z_'*X_) * Z_'*Z_ *inv(X_'*Z_);
stder_iv= sqrt(diag(varcoef)); % standerd errors
t_iv=beta_hat_iv./stder_iv; % t stats
//**************************************************************//
@#endif

@#if just_id==0
//***********IV estimation(general)*******************//
@#if nocon==0
x1_ = ones(sim_period__,1);
  //X_= [x1_ paipls_ mce_ ppai_];
 @#if ext_lg==1
 X_= [x1_  paipls_ mce_ ppai_ ppai2_ ppai3_ ppai4_];
 @#endif
 @#if ext_lg==0
 X_= [x1_  paipls_ mce_ ppai_ ];
 @#endif


 @#if mc_adj==0
 Z_= [x1_ fpaie_ mce_ ppai_];
 @#endif
 @#if mc_adj==1
 Z_= [x1_ ffpaie_ mce2_ ppai_];
 @#endif
 @#if mc_adj==2
  //Z_= [x1_ ffpaie_ fmce_ ppai_];
  @#if ext_lg==1
  Z_= [x1_ ffpaie_ fmce_ ppai_ ppai2_ ppai3_ ppai4_];
  @#endif
  @#if ext_lg==0
  Z_= [x1_ ffpaie_ fmce_ ppai_];
  @#endif
  
 @#endif
@#endif

@#if nocon==1
 @#if ext_lg==1
 X_= [ paipls_ mce_ ppai_ ppai2_ ppai3_ ppai4_];
 @#endif
 @#if ext_lg==0
 X_= [ paipls_ mce_ ppai_ ];
 @#endif


 @#if mc_adj==0
 Z_= [fpaie_ mce_ ppai_];
 @#endif
 @#if mc_adj==1
 Z_= [ffpaie_ mce2_ ppai_];
 @#endif
 @#if mc_adj==2
  @#if ext_lg==1
  Z_= [ffpaie_ fmce_ ppai_ ppai2_ ppai3_ ppai4_];
  @#endif
  @#if ext_lg==0
  Z_= [ffpaie_ fmce_ ppai_];
  @#endif
@#endif
@#endif

Pz_=Z_*(Z_'*Z_)^(-1)*Z_';
Xhat0_=Pz_*X_;

@#if nocon==0
  //Xhat_=Xhat0_(:,2:4);
  @#if ext_lg==1
  Xhat_=Xhat0_(:,2:7);
  @#endif
  @#if ext_lg==0
  Xhat_=Xhat0_(:,2:4);
  @#endif
@#endif

@#if nocon==1
  @#if ext_lg==1
  Xhat_=Xhat0_(:,1:6);
  @#endif
  @#if ext_lg==0
  Xhat_=Xhat0_(:,1:3);
  @#endif
@#endif

Y_ = pai_;
beta_hat_iv = ((X_'*Pz_*X_)^(-1))*X_'*Pz_*Y_; % IV estimated coeff
yhat = X_*beta_hat_iv;
res = Y_ - yhat; %residual ;
resvar = (res'*res)./(sim_period_);  %Var(residuals)
@#if nw_stdr==0
varcoef = resvar * inv(X_'*Pz_*X_);
stder_iv= sqrt(diag(varcoef)); % standerd errors
@#endif
@#if nw_stdr==1
 @#if nocon==0
 stder_iv=NeweyWest(res,Xhat_,L,1);
 @#endif
 @#if nocon==1
 stder_iv=NeweyWest(res,Xhat_,L,0);
 @#endif
@#endif

t_iv=beta_hat_iv./stder_iv; % t stats
//**************************************************************//
@#endif

@#endif


@#if merror==6
X=[fpai(1:sim_period_), mc(1:sim_period_), ppai(1:sim_period_)];
@#endif

@#if merror==5
X=[fpaie(1:sim_period_), mc(1:sim_period_), ppai(1:sim_period_)];
@#endif

@#if merror==4
X=[fpaie(1:sim_period_), mce(1:sim_period_), ppaie(1:sim_period_)];
@#endif

@#if merror==3
X=[fpaie(1:sim_period_), mce(1:sim_period_), ppai(1:sim_period_)];
@#endif

@#if merror==2
X=[fpaie(1:sim_period_), mce(1:sim_period_), ppai(1:sim_period_)];
@#endif

@#if merror==1
X=[fpai(1:sim_period_), mce(1:sim_period_), ppai(1:sim_period_)];
@#endif

@#if merror==0
X=[fpai(1:sim_period_), mc(1:sim_period_), ppai(1:sim_period_)];
@#endif

//@#if merror==7
//X=[pai(2:sim_period_+1), mce(1:sim_period_), ppai(1:sim_period_)];
//@#endif


//LinearModel.fit(X, pai_,'Intercept',false)
@#if merror==0
reg=LinearModel.fit(X, pai_);
@#endif
@#if merror==1
reg=LinearModel.fit(X, pai_);
@#endif
@#if merror==2
reg=LinearModel.fit(X, paiee_);
@#endif
@#if merror==3
reg=LinearModel.fit(X, pai_);
@#endif
@#if merror==4
reg=LinearModel.fit(X, paiee_);
@#endif
@#if merror==5
reg=LinearModel.fit(X, paiee_);
@#endif

@#if merror==6
reg=LinearModel.fit(X, paie_);
@#endif

@#if merror==7
//reg=LinearModel.fit(X_, tpai_);

@#if nocon==0
simseries_pccf{j,ii,1}=beta_hat_iv(1);
simseries_pccf{j,ii,2}=beta_hat_iv(2);
simseries_pccf{j,ii,3}=beta_hat_iv(3);
simseries_pccf{j,ii,4}=beta_hat_iv(4);

simseries_pcse{j,ii,1}=stder_iv(1);
simseries_pcse{j,ii,2}=stder_iv(2);
simseries_pcse{j,ii,3}=stder_iv(3);
simseries_pcse{j,ii,4}=stder_iv(4);

simseries_pct{j,ii,1}=t_iv(1);
simseries_pct{j,ii,2}=t_iv(2);
simseries_pct{j,ii,3}=t_iv(3);
simseries_pct{j,ii,4}=t_iv(4);
@#endif 

@#if nocon==1
simseries_pccf{j,ii,1}=beta_hat_iv(1);
simseries_pccf{j,ii,2}=beta_hat_iv(2);
simseries_pccf{j,ii,3}=beta_hat_iv(3);

simseries_pcse{j,ii,1}=stder_iv(1);
simseries_pcse{j,ii,2}=stder_iv(2);
simseries_pcse{j,ii,3}=stder_iv(3);

simseries_pct{j,ii,1}=t_iv(1);
simseries_pct{j,ii,2}=t_iv(2);
simseries_pct{j,ii,3}=t_iv(3);
@#endif 

@#endif

@#if merror==2
simseries_pccf{j,ii,1}=reg.Coefficients.Estimate(1);
simseries_pccf{j,ii,2}=reg.Coefficients.Estimate(2);
simseries_pccf{j,ii,3}=reg.Coefficients.Estimate(3);
simseries_pccf{j,ii,4}=reg.Coefficients.Estimate(4);

simseries_pcse{j,ii,1}=reg.Coefficients.SE(1);
simseries_pcse{j,ii,2}=reg.Coefficients.SE(2);
simseries_pcse{j,ii,3}=reg.Coefficients.SE(3);
simseries_pcse{j,ii,4}=reg.Coefficients.SE(4);

simseries_pct{j,ii,1}=reg.Coefficients.tStat(1);
simseries_pct{j,ii,2}=reg.Coefficients.tStat(2);
simseries_pct{j,ii,3}=reg.Coefficients.tStat(3);
simseries_pct{j,ii,4}=reg.Coefficients.tStat(4);

simseries_pcpv{j,ii,1}=reg.Coefficients.pValue(1);
simseries_pcpv{j,ii,2}=reg.Coefficients.pValue(2);
simseries_pcpv{j,ii,3}=reg.Coefficients.pValue(3);
simseries_pcpv{j,ii,4}=reg.Coefficients.pValue(4);
@#endif

//************************************************************************//
//** The End of PC curve estimates: 1=constant, 2=fpai, 3=gap(e), 4=ppai **//
//************************************************************************//

//*************************************************//
//** Variance, Std, Autocoeff(1st and 2nd lags) **//
//************************************************//



act_mcgap=autocorr(mcgap);
simseries_mcgap{j,ii,1}=var(mcgap);
simseries_mcgap{j,ii,2}=std(mcgap);
simseries_mcgap{j,ii,3}=act_mcgap(2);
simseries_mcgap{j,ii,4}=act_mcgap(3);

act_mc=autocorr(mc);
simseries_mc{j,ii,1}=var(mc);
simseries_mc{j,ii,2}=std(mc);
simseries_mc{j,ii,3}=act_mc(2);
simseries_mc{j,ii,4}=act_mc(3);

act_fmc=autocorr(fmc);
simseries_fmc{j,ii,1}=var(fmc);
simseries_fmc{j,ii,2}=std(fmc);
simseries_fmc{j,ii,3}=act_fmc(2);
simseries_fmc{j,ii,4}=act_fmc(3);

act_ffpai=autocorr(ffpai);
simseries_ffpai{j,ii,1}=var(ffpai);
simseries_ffpai{j,ii,2}=std(ffpai);
simseries_ffpai{j,ii,3}=act_ffpai(2);
simseries_ffpai{j,ii,4}=act_ffpai(3);

act_tfp=autocorr(tfp);
simseries_tfp{j,ii,1}=var(tfp);
simseries_tfp{j,ii,2}=std(tfp);
simseries_tfp{j,ii,3}=act_tfp(2);
simseries_tfp{j,ii,4}=act_tfp(3);

act_tfpu=autocorr(tfpu);
simseries_tfpu{j,ii,1}=var(tfpu);
simseries_tfpu{j,ii,2}=std(tfpu);
simseries_tfpu{j,ii,3}=act_tfpu(2);
simseries_tfpu{j,ii,4}=act_tfpu(3);

act_tfpn=autocorr(tfpn);
simseries_tfpn{j,ii,1}=var(tfpn);
simseries_tfpn{j,ii,2}=std(tfpn);
simseries_tfpn{j,ii,3}=act_tfpn(2);
simseries_tfpn{j,ii,4}=act_tfpn(3);

act_y=autocorr(y);
simseries_y{j,ii,1}=var(y);
simseries_y{j,ii,2}=std(y);
simseries_y{j,ii,3}=act_y(2);
simseries_y{j,ii,4}=act_y(3);

act_c=autocorr(c);
simseries_c{j,ii,1}=var(c);
simseries_c{j,ii,2}=std(c);
simseries_c{j,ii,3}=act_c(2);
simseries_c{j,ii,4}=act_c(3);

act_iv=autocorr(iv);
simseries_iv{j,ii,1}=var(iv);
simseries_iv{j,ii,2}=std(iv);
simseries_iv{j,ii,3}=act_iv(2);
simseries_iv{j,ii,4}=act_iv(3);

act_p=autocorr(p);
simseries_p{j,ii,1}=var(p);
simseries_p{j,ii,2}=std(p);
simseries_p{j,ii,3}=act_p(2);
simseries_p{j,ii,4}=act_p(3);

act_pai=autocorr(pai);
simseries_pai{j,ii,1}=var(pai);
simseries_pai{j,ii,2}=std(pai);
simseries_pai{j,ii,3}=act_pai(2);
simseries_pai{j,ii,4}=act_pai(3);

act_fpai=autocorr(fpai);
simseries_fpai{j,ii,1}=var(fpai);
simseries_fpai{j,ii,2}=std(fpai);
simseries_fpai{j,ii,3}=act_fpai(2);
simseries_fpai{j,ii,4}=act_fpai(3);

act_h=autocorr(h);
simseries_h{j,ii,1}=var(h);
simseries_h{j,ii,2}=std(h);
simseries_h{j,ii,3}=act_h(2);
simseries_h{j,ii,4}=act_h(3);

act_w=autocorr(w);
simseries_w{j,ii,1}=var(w);
simseries_w{j,ii,2}=std(w);
simseries_w{j,ii,3}=act_w(2);
simseries_w{j,ii,4}=act_w(3);

act_r=autocorr(r);
simseries_r{j,ii,1}=var(r);
simseries_r{j,ii,2}=std(r);
simseries_r{j,ii,3}=act_r(2);
simseries_r{j,ii,4}=act_r(3);

act_i=autocorr(i);
simseries_i{j,ii,1}=var(i);
simseries_i{j,ii,2}=std(i);
simseries_i{j,ii,3}=act_i(2);
simseries_i{j,ii,4}=act_i(3);

act_u=autocorr(u);
simseries_u{j,ii,1}=var(u);
simseries_u{j,ii,2}=std(u);
simseries_u{j,ii,3}=act_u(2);
simseries_u{j,ii,4}=act_u(3);


@#if tech==1

act_rd=autocorr(rd);
simseries_rd{j,ii,1}=var(rd);
simseries_rd{j,ii,2}=std(rd);
simseries_rd{j,ii,3}=act_rd(2);
simseries_rd{j,ii,4}=act_rd(3);

act_tiv=autocorr(tiv);
simseries_tiv{j,ii,1}=var(tiv);
simseries_tiv{j,ii,2}=std(tiv);
simseries_tiv{j,ii,3}=act_tiv(2);
simseries_tiv{j,ii,4}=act_tiv(3);

act_a=autocorr(a);
simseries_a{j,ii,1}=var(a);
simseries_a{j,ii,2}=std(a);
simseries_a{j,ii,3}=act_a(2);
simseries_a{j,ii,4}=act_a(3);

@#endif

//**********************************************************//
//** The End of Variance, Std, Autocoeff(1st and 2nd lags) **//
//**********************************************************//

//********************************************//
//** Correlation with fpai, fpai_t, fpai_dv***//
//********************************************//


cr_pai_fpaie=corrcoef(paipls_,ffpaie_);

@#if corr_spec_fpai==0
cr_fpai_fpaie=corrcoef(fpai,fpaie);
@#endif
@#if corr_spec_fpai==1
//cr_fpai_fpaie=corrcoef(ffpai(2:sim_period_+1),ffpaie_);
cr_fpai_fpaie=corrcoef(ffpai(cut_period1:sim_period_),ffpaie_);
@#endif

@#if corr_spec_mc==0
cr_mc_mce=corrcoef(mc,mce);
@#endif
@#if corr_spec_mc==1
cr_mc_mce=corrcoef(mc(cut_period2:sim_period_+1),mce_);
cr_mc_fmce=corrcoef(mc(cut_period2:sim_period_+1),fmce_);
cr_mce_fmce=corrcoef(mce_,fmce_);

@#endif


cr_fpaie_shk_meee=corrcoef(fpaie,gamma*shk_meee);
cr_fpaie_shk_me=corrcoef(fpaie,shk_me);
cr_fpaie_adjshk_meee=corrcoef(fpaie,(1-adj)*gamma*shk_meee);

cv_fpaie_shk_meee=cov(fpaie,gamma*shk_meee);
cv_fpaie_adjshk_meee=cov(fpaie,(1-adj)*gamma*shk_meee);

simseries_fpai_cr{j,ii,1}=cr_fpai_fpaie(2,1);
simseries_fpai_cr{j,ii,2}=cr_fpaie_shk_meee(2,1);
simseries_fpai_cr{j,ii,3}=cr_fpaie_shk_me(2,1);
simseries_fpai_cr{j,ii,4}=cr_fpaie_adjshk_meee(2,1);
simseries_fpai_cr{j,ii,5}=cr_pai_fpaie(2,1);

simseries_fpai_cv{j,ii,1}=cv_fpaie_shk_meee(2,1);
simseries_fpai_cv{j,ii,2}=cv_fpaie_adjshk_meee(2,1);

simseries_mc_mce_cr{j,ii,1}=cr_mc_mce(2,1);
simseries_mc_fmce_cr{j,ii,1}=cr_mc_fmce(2,1);
simseries_mce_fmce_cr{j,ii,1}=cr_mce_fmce(2,1);


//************************************//
//** End of Correlation with fpai***//
//***********************************//


//*****************************//
//** Correlation with Output**//
//** & Crr(ppai,pinf)  ***//
//****************************//

cr_y_c=corrcoef(y,c);
cr_y_iv=corrcoef(y,iv);
cr_y_h=corrcoef(y,h);
cr_y_p=corrcoef(y,p);
cr_y_pai=corrcoef(y,pai);
cr_y_i=corrcoef(y,i);
cr_y_r=corrcoef(y,r);
cr_y_w=corrcoef(y,w);
cr_y_u=corrcoef(y,u);
cr_y_tfp=corrcoef(y,tfp);

@#if tech==1
cr_y_a=corrcoef(y,a);
cr_y_rd=corrcoef(y,rd);
cr_y_tiv=corrcoef(y,tiv);
cr_ppai_pinf=corrcoef(ppai,pinf);
@#endif

simseries_cr{j,ii,1}=cr_y_c(2,1);
simseries_cr{j,ii,2}=cr_y_iv(2,1);
simseries_cr{j,ii,3}=cr_y_h(2,1);
simseries_cr{j,ii,4}=cr_y_p(2,1);
simseries_cr{j,ii,5}=cr_y_pai(2,1);
simseries_cr{j,ii,6}=cr_y_i(2,1);
simseries_cr{j,ii,7}=cr_y_r(2,1);
simseries_cr{j,ii,8}=cr_y_w(2,1);
simseries_cr{j,ii,9}=cr_y_u(2,1);
simseries_cr{j,ii,10}=cr_y_tfp(2,1);

@#if tech==1
simseries_cr{j,ii,11}=cr_y_a(2,1);
simseries_cr{j,ii,12}=cr_y_rd(2,1);
simseries_cr{j,ii,13}=cr_y_tiv(2,1);
simseries_cr{j,ii,14}=cr_ppai_pinf(2,1);
@#endif

//************************************//
//** End of Correlation with Output***//
//***********************************//

//**end of fig==0**//
@#endif
//****************//


end
//** End of 'for ii' **//
end 
//** End of 'for j' **//



//** start of fig==0 **//
@#if fig==0  
//********************//

//*****************************//
//** Var, STD and  AutoCorr***//
//****************************//
std_mcgap=[simseries_mcgap{:,:,2}];
std_mc=[simseries_mc{:,:,2}];
std_fmc=[simseries_fmc{:,:,2}];

var_tfp=[simseries_tfp{:,:,1}];
std_tfp=[simseries_tfp{:,:,2}];
acrr1_tfp=[simseries_tfp{:,:,3}];
acrr2_tfp=[simseries_tfp{:,:,4}];

var_tfpu=[simseries_tfpu{:,:,1}];
std_tfpu=[simseries_tfpu{:,:,2}];
acrr1_tfpu=[simseries_tfpu{:,:,3}];
acrr2_tfpu=[simseries_tfpu{:,:,4}];

var_tfpn=[simseries_tfpn{:,:,1}];
std_tfpn=[simseries_tfpn{:,:,2}];
acrr1_tfpn=[simseries_tfpn{:,:,3}];
acrr2_tfpn=[simseries_tfpn{:,:,4}];


var_y=[simseries_y{:,:,1}];
std_y=[simseries_y{:,:,2}];
acrr1_y=[simseries_y{:,:,3}];
acrr2_y=[simseries_y{:,:,4}];

var_c=[simseries_c{:,:,1}];
std_c=[simseries_c{:,:,2}];
acrr1_c=[simseries_c{:,:,3}];
acrr2_c=[simseries_c{:,:,4}];

var_iv=[simseries_iv{:,:,1}];
std_iv=[simseries_iv{:,:,2}];
acrr1_iv=[simseries_iv{:,:,3}];
acrr2_iv=[simseries_iv{:,:,4}];

var_h=[simseries_h{:,:,1}];
std_h=[simseries_h{:,:,2}];
acrr1_h=[simseries_h{:,:,3}];
acrr2_h=[simseries_h{:,:,4}];

var_p=[simseries_p{:,:,1}];
std_p=[simseries_p{:,:,2}];
acrr1_p=[simseries_p{:,:,3}];
acrr2_p=[simseries_p{:,:,4}];

var_pai=[simseries_pai{:,:,1}];
std_pai=[simseries_pai{:,:,2}];
acrr1_pai=[simseries_pai{:,:,3}];
acrr2_pai=[simseries_pai{:,:,4}];

var_fpai=[simseries_fpai{:,:,1}];
std_fpai=[simseries_fpai{:,:,2}];
acrr1_fpai=[simseries_fpai{:,:,3}];
acrr2_fpai=[simseries_fpai{:,:,4}];

var_ffpai=[simseries_ffpai{:,:,1}];
std_ffpai=[simseries_ffpai{:,:,2}];
acrr1_ffpai=[simseries_ffpai{:,:,3}];
acrr2_ffpai=[simseries_ffpai{:,:,4}];

var_i=[simseries_i{:,:,1}];
std_i=[simseries_i{:,:,2}];
acrr1_i=[simseries_i{:,:,3}];
acrr2_i=[simseries_i{:,:,4}];

var_r=[simseries_r{:,:,1}];
std_r=[simseries_r{:,:,2}];
acrr1_r=[simseries_r{:,:,3}];
acrr2_r=[simseries_r{:,:,4}];

var_w=[simseries_w{:,:,1}];
std_w=[simseries_w{:,:,2}];
acrr1_w=[simseries_w{:,:,3}];
acrr2_w=[simseries_w{:,:,4}];

var_u=[simseries_u{:,:,1}];
std_u=[simseries_u{:,:,2}];
acrr1_u=[simseries_u{:,:,3}];
acrr2_u=[simseries_u{:,:,4}];

@#if tech==1
var_a=[simseries_a{:,:,1}];
std_a=[simseries_a{:,:,2}];
acrr1_a=[simseries_a{:,:,3}];
acrr2_a=[simseries_a{:,:,4}];

var_rd=[simseries_rd{:,:,1}];
std_rd=[simseries_rd{:,:,2}];
acrr1_rd=[simseries_rd{:,:,3}];
acrr2_rd=[simseries_rd{:,:,4}];

var_tiv=[simseries_tiv{:,:,1}];
std_tiv=[simseries_tiv{:,:,2}];
acrr1_tiv=[simseries_tiv{:,:,3}];
acrr2_tiv=[simseries_tiv{:,:,4}];
@#endif


//*********************************************//
//** Correlation with Output &crr(ppai,pinf)***//
//*******************************************//
crr_y_c=[simseries_cr{:,:,1}];
crr_y_iv=[simseries_cr{:,:,2}];
crr_y_h=[simseries_cr{:,:,3}];
crr_y_p=[simseries_cr{:,:,4}];
crr_y_pai=[simseries_cr{:,:,5}];
crr_y_i=[simseries_cr{:,:,6}];
crr_y_r=[simseries_cr{:,:,7}];
crr_y_w=[simseries_cr{:,:,8}];
crr_y_u=[simseries_cr{:,:,9}];
crr_y_tfp=[simseries_cr{:,:,10}];

@#if tech==1
crr_y_a=[simseries_cr{:,:,11}];
crr_y_rd=[simseries_cr{:,:,12}];
crr_y_tiv=[simseries_cr{:,:,13}];
crr_ppai_pinf=[simseries_cr{:,:,14}];
@#endif




//*******************//
//** PC estimates***//
//******************//

//*** With constant terms ***//
@#if nocon==0
pccf2=[simseries_pccf{:,:,2}];
pccf3=[simseries_pccf{:,:,3}];
pccf4=[simseries_pccf{:,:,4}];
@#endif
//*** NO constant terms ***//
@#if nocon==1
pccf2=[simseries_pccf{:,:,1}];
pccf3=[simseries_pccf{:,:,2}];
pccf4=[simseries_pccf{:,:,3}];
@#endif


//*** the following command produces the index of corrsponding value(e.g., 'im_pccf2' gives the index which has the median value in pccf2 vector)***//  
[cm_pccf2,im_pccf2]=min((pccf2-median(pccf2)).^2);
[cm_pccf3,im_pccf3]=min((pccf3-median(pccf3)).^2);
[cm_pccf4,im_pccf4]=min((pccf4-median(pccf4)).^2);


[c1_pccf2,i1_pccf2]=min((pccf2-prctile(pccf2,1)).^2);
[c1_pccf3,i1_pccf3]=min((pccf3-prctile(pccf3,1)).^2);
[c1_pccf4,i1_pccf4]=min((pccf4-prctile(pccf4,1)).^2);


[c2p5_pccf2,i2p5_pccf2]=min((pccf2-prctile(pccf2,2.5)).^2);
[c2p5_pccf3,i2p5_pccf3]=min((pccf3-prctile(pccf3,2.5)).^2);
[c2p5_pccf4,i2p5_pccf4]=min((pccf4-prctile(pccf4,2.5)).^2);

[c5_pccf2,i5_pccf2]=min((pccf2-prctile(pccf2,5)).^2);
[c5_pccf3,i5_pccf3]=min((pccf3-prctile(pccf3,5)).^2);
[c5_pccf4,i5_pccf4]=min((pccf4-prctile(pccf4,5)).^2);

[c10_pccf2,i10_pccf2]=min((pccf2-prctile(pccf2,10)).^2);
[c10_pccf3,i10_pccf3]=min((pccf3-prctile(pccf3,10)).^2);
[c10_pccf4,i10_pccf4]=min((pccf4-prctile(pccf4,10)).^2);

[c20_pccf2,i20_pccf2]=min((pccf2-prctile(pccf2,20)).^2);
[c20_pccf3,i20_pccf3]=min((pccf3-prctile(pccf3,20)).^2);
[c20_pccf4,i20_pccf4]=min((pccf4-prctile(pccf4,20)).^2);


[c30_pccf2,i30_pccf2]=min((pccf2-prctile(pccf2,30)).^2);
[c30_pccf3,i30_pccf3]=min((pccf3-prctile(pccf3,30)).^2);
[c30_pccf4,i30_pccf4]=min((pccf4-prctile(pccf4,30)).^2);

[c40_pccf2,i40_pccf2]=min((pccf2-prctile(pccf2,40)).^2);
[c40_pccf3,i40_pccf3]=min((pccf3-prctile(pccf3,40)).^2);
[c40_pccf4,i40_pccf4]=min((pccf4-prctile(pccf4,40)).^2);

[c50_pccf2,i50_pccf2]=min((pccf2-prctile(pccf2,50)).^2);
[c50_pccf3,i50_pccf3]=min((pccf3-prctile(pccf3,50)).^2);
[c50_pccf4,i50_pccf4]=min((pccf4-prctile(pccf4,50)).^2);

[c60_pccf2,i60_pccf2]=min((pccf2-prctile(pccf2,60)).^2);
[c60_pccf3,i60_pccf3]=min((pccf3-prctile(pccf3,60)).^2);
[c60_pccf4,i60_pccf4]=min((pccf4-prctile(pccf4,60)).^2);

[c70_pccf2,i70_pccf2]=min((pccf2-prctile(pccf2,70)).^2);
[c70_pccf3,i70_pccf3]=min((pccf3-prctile(pccf3,70)).^2);
[c70_pccf4,i70_pccf4]=min((pccf4-prctile(pccf4,70)).^2);

[c80_pccf2,i80_pccf2]=min((pccf2-prctile(pccf2,80)).^2);
[c80_pccf3,i80_pccf3]=min((pccf3-prctile(pccf3,80)).^2);
[c80_pccf4,i80_pccf4]=min((pccf4-prctile(pccf4,80)).^2);

[c90_pccf2,i90_pccf2]=min((pccf2-prctile(pccf2,90)).^2);
[c90_pccf3,i90_pccf3]=min((pccf3-prctile(pccf3,90)).^2);
[c90_pccf4,i90_pccf4]=min((pccf4-prctile(pccf4,90)).^2);

[c95_pccf2,i95_pccf2]=min((pccf2-prctile(pccf2,95)).^2);
[c95_pccf3,i95_pccf3]=min((pccf3-prctile(pccf3,95)).^2);
[c95_pccf4,i95_pccf4]=min((pccf4-prctile(pccf4,95)).^2);

[c97p5_pccf2,i97p5_pccf2]=min((pccf2-prctile(pccf2,97.5)).^2);
[c97p5_pccf3,i97p5_pccf3]=min((pccf3-prctile(pccf3,97.5)).^2);
[c97p5_pccf4,i97p5_pccf4]=min((pccf4-prctile(pccf4,97.5)).^2);


[c99_pccf2,i99_pccf2]=min((pccf2-prctile(pccf2,99)).^2);
[c99_pccf3,i99_pccf3]=min((pccf3-prctile(pccf3,99)).^2);
[c99_pccf4,i99_pccf4]=min((pccf4-prctile(pccf4,99)).^2);



//**************************************************************//

@#if nocon==0
pct2=[simseries_pct{:,:,2}];
pct3=[simseries_pct{:,:,3}];
pct4=[simseries_pct{:,:,4}];
@#endif
@#if nocon==1
pct2=[simseries_pct{:,:,1}];
pct3=[simseries_pct{:,:,2}];
pct4=[simseries_pct{:,:,3}];
@#endif

@#if nocon==0
pcse2=[simseries_pcse{:,:,2}];
pcse3=[simseries_pcse{:,:,3}];
pcse4=[simseries_pcse{:,:,4}];
@#endif
@#if nocon==1
pcse2=[simseries_pcse{:,:,1}];
pcse3=[simseries_pcse{:,:,2}];
pcse4=[simseries_pcse{:,:,3}];
@#endif



crr_fpai_fpaie=[simseries_fpai_cr{:,:,1}];
crr_fpaie_shk_meee=[simseries_fpai_cr{:,:,2}];
crr_fpaie_shk_me=[simseries_fpai_cr{:,:,3}];
crr_fpaie_adjshk_meee=[simseries_fpai_cr{:,:,4}];
crr_pai_fpaie=[simseries_fpai_cr{:,:,5}];

crr_mc_mce=[simseries_mc_mce_cr{:,:,1}];
crr_mc_fmce=[simseries_mc_fmce_cr{:,:,1}];
crr_mce_fmce=[simseries_mce_fmce_cr{:,:,1}];

cov_fpaie_shk_meee=[simseries_fpai_cv{:,:,1}];
cov_fpaie_adjshk_meee=[simseries_fpai_cv{:,:,2}];

corr_mc1_nm=['corr(mc, mce)']
corr_mc_mce=[mean(crr_mc_mce)]


pccf4_max= pccf4(i99_pccf4);
pccf4_min=pccf4(i1_pccf4);
pccf2_max= pccf2(i99_pccf2);
pccf2_min=pccf2(i1_pccf2);


mtr=[pccf4' pccf2' pccf3' pct4'  pct2'  pct3'  pcse4'  pcse2'  pcse3'];
dat1=mtr;
idx_pc2max=dat1(:, 2)>=pccf2_max;
dat2=dat1(~idx_pc2max,:);% excluding rows which has >=pccf_max
idx_pc2min=dat2(:, 2)<=pccf2_min;
dat3=dat2(~idx_pc2min,:);% excluding rows which has <=pccf_min
idx_pc4max=dat3(:, 1)>=pccf4_max;
dat4=dat3(~idx_pc4max,:);% excluding rows which has >=pccf_max
idx_pc4min=dat4(:, 1)<=pccf4_min;
dat=dat4(~idx_pc4min,:);% excluding rows which has <=pccf_min




//** end of fig==0 **//
@#endif
//**********************//

